All-electron magnetic response with pseudopotentials: NMR 

chemical shifts 



Chris J. PickardB 

Institut fur Geowissenschaften, 
Universitdt Kiel, Olshausenstrasse 40 D-24-098 Kiel, Germany 

Francesco Mauri 
Laboratoire de Mineralogie-Cristallographie de Paris, 
Universite Pierre et Marie Curie, 4 Place Jussieu, 75252, Paris, Cedex 05, France 

(February 1, 2008) 

Abstract 

A theory for the ab initio calculation of all-electron NMR chemical shifts 
in insulators using pseudopotentials is presented. It is formulated for both 
finite and infinitely periodic systems and is based on an extension to the 
Projector Augmented Wave approach of Blochl [P. E. Blochl, Phys. Rev. B 
50, 17953 (1994)] and the method of Mauri et al [F. Mauri, B. G. Pfrommer, 
and S. G. Louie, Phys. Rev. Lett. 77, 5300 (1996)]. The theory is successfully 
validated for molecules by comparison with a selection of quantum chemical 
results, and in periodic systems by comparison with plane-wave all-electron 
results for diamond. 
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I. INTRODUCTION 



The experimental technique of nuclear magnetic resonance (NMR) is widely used in 
structural chemistry and increasingly in solid state studies.0 Chemical shift (a) spectra give 
information about the atomic structure of the sample under investigation. In the case of 
molecular systems, empirical rules are commonly used to extract this information from 
the raw experimental data. However, this approach cannot be applied in the solid state, 
as the atomic configurations often cannot be modeled by chemical analogues or reference 
compounds. In these cases ab initio calculations of the chemical shifts are the only way to 
obtain an unambiguous determination of the microscopic structure. 

Until recently, there has been no theory for the calculation of NMR chemical shifts in 
extended periodic systems, and the conventional approach to the theoretical interpretation 
of solid state NMR spectra has been to approximate the infinite solid by a cluster.! In this 
way, the traditional quantum chemical approaches! - ! can be used to calculate the chemical 
shifts. Unfortunately, true convergence with respect to basis set and cluster size is often not 
possible due to the limitations of available computational resources. 

The work of Mauri, Pfrommer and Louie! solved the problem of calculating NMR chemi- 
cal shifts in the solid state with an all-electron Hamiltonian. Integrated with their approach 
to the calculation of magnetic susceptibility,! they presented a theory for the ab initio com- 
putation of NMR chemical shifts in condensed matter systems using periodic boundary 
conditions (hereafter referred to as the MPL method). Although the MPL theory has been 
derived using an all-electron Hamiltonian, so far, it has only been implemented in an elec- 
tronic structure code based on norm-conserving pseudopotentials. In such implementation 
the complications inherent within the pseudopotential approximation have been neglected. 
For this reason, while several useful applications have emerged,IHi3 the method's use has 
been restricted to the calculation of chemical shifts of light elements (hydrogen, carbon, and 
nitrogen) and of silicon. Moreover, the description of the silicon chemical shifts required 
the explicit inclusion of the 2s and 2p silicon orbitals as valence and the use of a very high, 
and computationally expensive, plane-wave energy- cutoff of 600 Ry.0 In the above appli- 
cations of the MPL method the pseudopotential error had been assumed to be small and 
controllable. To compute the NMR chemical shifts of nuclei heavier than neon and to truly 
exploit the ability of pseudopotentials to calculate the properties of complex, low symmetry 
structures (which is well established for a wide range of structural properties), a theory is 
required which does not ignore the pseudopotential approximation. 

Apart from the early and isolated attempt of Ridard, Levy and Milliejjl it has been 
widely expected within the quantum chemical community that any theory for the calculation 
of NMR chemical shifts for nuclei described with a pseudopotential^ would fail due to the 
non-rigid nature of the core contributions to the total chemical shift.! However, a careful 
separation of core and valence contributions that ensures that they are individually gauge- 
invariant, by Gregor, Mauri and CarJlZl has shown that this is not the case and that the core 
contributions are rigid. This suggests that a pseudopotential based theory of NMR might, 
in fact, exist. 

One of the most obvious deficiencies of the pseudopotential approach is that the pseu- 
dopotential approximation explicitly neglects the form of the electronic wavefunctions near 
the nucleus. The pseudo wavefunctions are chosen to be as smooth as possible in the core 
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region, and the correct nodal structure of the wavefunctions is lost. This leads to a good 
approximation for the calculation of total energies and their derivatives, and properties 
for which the matrix elements are dominated by the regions outside the core. However, 
the quantitative calculation of many properties — hyperfine parameters, core level spectra, 
electric-field-gradients and the NMR chemical shifts — depend critically on the details of 
the all-electron wavefunctions at the nucleus. Van de Walle and Blochl presented a solution 
to this problem for the calculation of hyperfine parameters^ based on Blochl's Projector 
Augmented- Wave (PAW) electronic structure methodjl! which is itself closely related to 
Vanderbilt's ultrasoft pseudopotential schemeS While in all but a few reported cases, where 
core-electron polarization effects are importantlll or in some magnetic systems ,0 the PAW 
method gives similar results to pseudopotential approaches, it does provide a extremely 
useful framework for the unification of all-electron (Full-potential) Linearized Augmented 
Plane- Waveil and pseudopotential approaches. Indeed, it it becoming clear that the PAW 
approach, which will be described in more detail in Section 111 A , offers a general approach 
to the calculation of all-electron properties from pseudopotential based schemes. Following 
the work of Van de Walle and Blochl, core level spectra,il'i! momentum matrix elementsjl 
and electric field gradients^ have all been calculated using the PAW scheme. 

In this paper we present a theory for all-electron magnetic response within the pseudopo- 
tential approximation and its application to the calculation of first principles NMR chemical 
shifts. The connection between the current response and the chemical shifts is outlined in 
Section ||. We introduce an extension of Blochl's PAW approach, which we call the Gauge 
Including Projector Augmented- Wave (GIPAW) approach. This will be described in Section 
III. A Hamiltonian constructed using GIPAW has the required translational invariance in 
the presence of a magnetic field. This is not true for the original PAW formulation. In 
Section [TV] we present our theory for finite systems. In Section [V] we reformulate our expres- 
sions for extended systems. To be useful, these expressions must be restricted to periodic 
extended systems, and the periodic theory is presented in Section |V|. Both the theories for 
finite and extended periodic systems summarized in Section [VII] have been implemented in 
a plane-wave pseudopotential electronic structure code. Details of our implementation are 
given in Section |V1U| . We validate the method by comparison with calculations by Gregor et 
aJll for a selection of small molecules. The theory for extended systems is further validated 
by comparison to results obtained by an all-electron plane-wave calculation for a crystalline 
material, diamond. 



II. NMR CHEMICAL SHIFTS 

A uniform, external magnetic field B applied to a sample of matter induces an electric 
current. In an insulating non-magnetic material, only the orbital motion of the electrons 
contribute to this current. Moreover, for the field strengths typically used in NMR experi- 
ments, the induced electronic current is proportional to the external field B and is the first 
order induced current, j^(r'). The current j^(r') produces a non- uniform magnetic field, 

B«(r) = l/rfVjW(r')x^-. (1) 
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The chemical shift is defined as the ratio between the induced magnetic field and the external 
uniform applied magnetic field: 

B«(r) = -a(r)B. (2) 

Here c?(r) is the chemical shift tensor, and the isotropic chemical shift is given by a(r) = 
Tr[(j(r)]/3. NMR experiments can measure <r(r) at the nuclear positions. To compute the 
chemical shift tensor we first obtain j^(r) by perturbation theory and then we evaluate 
B-p (r) using Eq. ([[]). We now describe our new approach to the calculation of an induced 
all-electron current j^(r') using pseudopotentials and in Section |V11J| the computational 
procedure we use to obtain j^r'), and finally cr(r), is detailed. 



III. PSEUDOPOTENTIALS IN A MAGNETIC FIELD 

In this section we develop the Gauge Including Plane- Wave method, first describing the 
original Projector Augmented- Wave method, and then extending it to the case of a uniform 
applied magnetic field. 

A. Projector augmented-wave method 



In Ref. [L^, Blochl introduced a linear transformation operator T which maps the valence 
pseudo wavefunctions |\&) onto the corresponding all-electron wavef unctions, = T\^f). 
The operator is defined by specifying a set of target all-electron partial waves \<f>R in ) obtained 
by the application of T on to a set of pseudo partial waves |0R, n ) with 

T=l + ^[i0 Rin )-|0 Rin )](pK,n| (3) 
R,n 

and (p-R t n\ are a set of projectors such that (pR, n |0R',m) — ^R,R'^n,m- Each projector and 
partial wave is an atomic-like function centered on an atomic site R, and the index n refers 
to the angular momentum quantum numbers and to an additional number used if there 
are more than one projector per angular momentum channel. The expectation value of an 
operator O between all-electron wavefunctions can be expressed as the expectation value of 
a pseudo operator O = T + OT between the corresponding pseudo wavefunctions. 

To obtain a useful formalism we must make some further assumptions. In particular, 
for each atomic site we define an augmentation region f2 R and suppose that: i) outside the 
augmentation region the |0R )fl ) coincide with the |0r )TI ), ii) outside the augmentation 
region Q R , the |p R n ) vanish, iii) within the augmentation region fi R , the |0R, n ) form a com- 
plete set for the valence wavefunctions, i.e. any physical valence all-electron wavefunction 
can be written, within f2 R , as a linear combination of all-electron partial waves, and finally 
iv) the augmentation regions of different sites do not overlap. Blochl has shown that given 
these assumptions, if O is a local or a semi-local (such as p or p 2 ) operator: 

= 0+ |PR,n)[(0R,n|O |0R,m) — (0R,n \0\4>B.,m)}(pK,m\- (4) 

R,n,m 
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For simplicity, we shall further suppose that the norms computed within flu of \4>n,n) 
and |0R, n ) coincide. We then recover the norm conserving pseudopotential formalism in 
the Kleinman-Bylanderl! form. The pseudo wavefunctions which correspond to the all- 
electron valence eigenstates of the all-electron Hamiltonian H are eigenstates of the pseudo 
Hamiltonian H with the same eigenvalues. In the absence of a magnetic field the pseudo 
Hamiltonian is: 

H = T+HT = ip 2 + V loc (r) + £ Vg, (5) 

1 R 

where p is the momentum operator, and \^ loc (r) is the local part of the pseudopotentials, 
which includes the self-consistent part of the Hamiltonian. The non-local part of the pseu- 
dopotential at the atomic site R in the above expression is, 

t# = El*W<m<PR.J- (6) 

n,m 

The a^ m are the strengths of the non-local potential in each channel, and the depend on R 
since each atomic site may be occupied by a different chemical species. 

The choice of the pseudo partial waves and projectors is largely arbitrary. However, for a 
scheme to be useful, all the lowest eigenvalues of H should coincide with a valence eigenvalue 
of H up to an given energy E™ x , i.e. no ghost states should be introduced in the pseudo 
spectrum up to an energy -E™f x . The energy -E™f x depends on the specific property we wish 
to compute, and should at least be larger than the highest occupied eigenvalue. 

In contrast to the traditional formulation of pseudopotentials, using the PAW formulation 
it is possible to obtain the expectation values of all-electron operators in terms of pseudo 
wavefunctions using the pseudo operators defined in Eq. (f|). 



B. A single augmentation region in a uniform magnetic field 

In presence of a uniform external magnetic field B the all-electron Hamiltonian is: 

H = \(p + - c Mr)) 2 + V(r), (7) 

where c is the speed of light, V^(r) is the all-electron local potential, and B = V x A(r). 
We want to construct the corresponding pseudo Hamiltonian for a complex system, which 
will contain many augmentation regions. However, before treating this general case, we 
consider a simplified system with just a single augmentation region. The spatial origin is 
chosen to coincide with the atomic site of the augmentation region. In the symmetric gauge 
A(r) = ~B x (r — d), where d is a constant vector which indicates the gauge origin. The 
expectation values of the all-electron eigenstates for observable operators do not depend on 
the gauge origin d. However, the number of partial waves required to correctly describe the 
valence all-electron eigenstates in the augmentation region critically depends on the choice 
of d. To minimize the number of partial waves required we must put the gauge origin at 
the atomic site of the augmentation region, setting d = 0. Making this choice, we minimize 
the effect of the magnetic field on the all-electron wavefunctions in the augmentation region, 
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where |A(r)| 2 and its spatial derivatives attain their minimum value. Moreover, with this 
choice of gauge, the interaction between the valence and core states of the augmented atom 
is negligibly small.0 This is essential if we are to make the pseudopotential approximation. 
With 

A(r) = x r, (8) 

the all-electron Hamiltonian becomes: 

H = \v 2 + V(v) + II • B + ^(B x if, (9) 

where L = r x p is the angular momentum operator computed with respect to the atomic site 
within the augmentation region. Using Eq. (|J) and Eq. we obtain the corresponding 
pseudo Hamiltonian: 

H = \v 2 + V loc (r) + V» l + i-L • B + -L(B x if + £ \p , n )(b^ + b^)(p , m \, (10) 

Z AC oC nm 

where 

b ( nl = • [(<ML|0O, m > - (4> ,n\M4> ,m)} (n) 

and 

&£L = A[(^o,n|(B x if |0 o , m > - (0o,n|(B x if |0 o , m >] (12) 

If just one projector per angular momentum channel is used, as is usually the case with norm 
conserving pseudopotentialsJiH b^ m exactly vanishes, since |0o,n) an d |0o,n) are eigenstates 
of L and with the same norm within the augmentation region. Moreover, since (B x rf 
goes to zero in the center of the augmentation region, for norm conserving pseudopoten- 
tials the term b^ m can also be neglected. Thus, with one augmentation region centered 
at the gauge origin, the coupling with the magnetic field in the pseudo and all-electron 
Hamiltonians has the same form, i.e.: 

H = \v 2 + V lo \v) + V? + ±L ■ B + ^(B x if. (13) 



C. Translations in a uniform magnetic field 

The derivation in the previous section is not useful for systems with several augmentation 
regions. Indeed, the gauge origin can coincide with just one augmentation site at any given 
time. As a result, the number for projectors of the other augmentation regions would have 
to be increased to reach completeness in those regions. The cause of this problem is that 
the PAW approach does not preserve translational invariance in a uniform magnetic field. 

In a uniform magnetic field the description of the system should be invariant upon a 
rigid translation of all the atoms by a vector t. Following the translation, the all-electron 
potential becomes V'(r) = V(r — t) and the corresponding Hamiltonian is: 
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iT=i(p + ^A(r)) 2 + y(r-t), (14) 

where A(r) is still given by Eq. (|8|). Because of the translational invariance, the eigenenergies 
of H' coincide with the eigenenergies of the original Hamiltonian H. However, the new 
eigenstates l^), are not just obtained by a rigid translation of the original eigenstates \^ n ), 
but, upon translation, they pick up an additional phase factor proportional to the magnetic 
field: 

(r|^) = e^ txB (r-t|vI/ n ). (15) 

The PAW transformation does not ensure exact invariance upon translation, since the pseudo 
wavefunctions constructed with the T transformation operator of Eq. (Q) do not transform 
according to Eq. (|l5l) . 



D. Gauge including projector augmented-wave method 

To restore the translational invariance within a PAW-like approach, we introduce a field 
dependent transformation operator 7b, which, by construction, imposes the translational 
invariance exactly: 

T B = 1 + £^ r - RXB [|0R,n> _ |^ R)n)]( £ R)n | e -£r.RxB_ 
R,n 

This new transformation defines our novel approach, which we call the Gauge Including 
Projected Augmented- Wave (GIPAW) method. In the following, we indicate with a bar the 
pseudo wavefunctions and operators obtained using 7b operator by analogy to Blochl's use 
of the tilde. By construction, the pseudo eigenstates, generated from the all-electron 
eigenstates using \^/) = 7b \^f), satisfy the same translation relation as the all-electron 
eigenstates given by Eq. fllSf) . The GIPAW pseudo operator O = T^OT B corresponding to 
a local or a semi-local operator O is given by: 

= 0+ £ e^ RxB |p R , n > [(0R,n|e-^ r - RxB Oe^ r - RxB |0 R , m ) - (0R,,,|e-^ r - RxB Oe^ r - RxB |0 R , m ); 

R,n,m 

(PR, m |e-^ rRxB (17) 

There are connections between our GIPAW approach, and the gauge-including atomic 
orbitalsEI (GIAO) and the independent gauge for localized orbitalsi (IGLO) methods, widely 
used in the quantum chemical community. However, it should be recognized that in GI- 
PAW the phase required to maintain the translational invariance is carried by the operators, 
whereas in the GIAO and in IGLO approaches the field dependent phase is attached to the 
basis functions and to the occupied electronic orbitals, respectively. 



E. GIPAW Hamiltonian 

Using Eq. (fi7l), the identity 
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e -£r.RxB 



(p + ^A(r))" e^ rRxB = (p + ^A(r - R))™ , (18) 
for integer n, and the outcomes of the discussion concerning and 6^ in Section |III B 



we finally obtain the GIPAW pseudo Hamiltonian: 

H = -p 2 + V loc (r) + £ e^ r - RxB Vg 1 e-* r - RxB + — L ■ B + -^(B x r) 2 . (19) 

The GIPAW Hamiltonian coincides with the PAW Hamiltonian, Eq. (|j) for B = 0, and 
with the PAW Hamiltonian, Eq. (|T3|), for B ^ in systems with a single augmentation 
region centered at the origin. Moreover, as expected, the GIPAW eigenenergies are exactly 
invariant upon translation, in contrast to the PAW eigenenergies. 

For later use in perturbation theory, H can be expanded in powers of B\ 

H = # (0) + + 0(B 2 ) (20) 
where = is the unperturbed Hamiltonian given by Eq. (El), and 



^=^(L + i;Rxvg).B, 



(21) 



where 

1 



R- 7 [r,<], (22) 

t 



and with square brackets we indicate the commutator. 



F. GIPAW Current operator 

Another observable required to compute the NMR chemical shifts is the current. The 
all-electron electric current operator evaluated at the position r' is: 

J(r') = JV) - ^|r')(r'| = .P(r') - ^V>(r'|, (23) 
where J p (r') is the paramagnetic current operator, 

JV) = - P|r/)(r/ ' + |r/)(r/|P . (24) 
Using Eq. (|17D, and Eq. flT8|), we obtain the corresponding GIPAW operator: 



J(r') = J p (r') - ^^V) (r'| + ]T e £ rRxB [AJ^(r') + AJ^(r')] e -^ r - RxB , (25) 

2C R 

where 

AJ R (r') = ]T |pa, n ) \{W, n \^{v')\4>n,m) ~ (0R,n|J p (r')|0R, m )l (Pn, m \ (26) 
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is what we call the paramagnetic augmentation operator, and 
B x (V - R) 



AJ R (r') 



2c 



13 |PR,n) (0R,n|r')(r'|0R,m) - (0R,n|r')(r'|^R, m ) (pR, m |, (27) 



is what we call the diamagnetic augmentation operator. 

As for the Hamiltonian, for perturbation theory purposes it is useful to expand the 
operator J(r) in powers of B: 



with 



and 



jW(r') 



J(r') = J (0) (r') + 3 {1 \r') + 0(B< 



j( )(r') = J p (rO + EAJR(r') 

R 



B X r '\r')(v'\ + £ [AJ* (r') + ^[B x R ■ r, AJ^(r')] 



2c 



(28) 



(29) 



(30) 



IV. CURRENT RESPONSE IN FINITE SYSTEMS 

Within density functional perturbation theory, the current can be computed using the 
GIPAW operators and wavefunctions as: 

jW(r') = 2J2 [<*?>| J<°>(OW 0) > + (W (0) C)I*? ) > + (^ 0) |J (1) (r')l^ 0) >] • (31) 

o 

Here the factor of two accounts for spin degeneracy and the sum runs over the occupied 
orbitals o. The wavefunction \^$') is an unperturbed eigenstate of with eigenvalue e n 
and \^f^) is its linear variation, projected in the empty subspace: 

\^ ] ) = G(s n )H^\^). (32) 

The Green function operator is: 

m=E \mm, (33 ) 

with the sum running over the empty orbitals e. Reordering the different contributions of 
Eq. fl3"ip we obtain: 

j (1 HrO=jS rc (r')+jS(r / )+jk i ; i (r / ), (34) 

where 

j£L(0 = 4 H Re [(^inr'Meo)^^)} - ^P PS (r')B X r'. (35) 
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Re stands for taking the real part and p ps (r') = 2 J2o(^^\ r> ) ( v '\^ ( o' > ) * s the ground state 
pseudo density. The paramagnetic correction to the current is 



R' o 



and the diamagnetic correction is 

£ 1 kr') = 2E(^° ) |AJ d R(r')l^ 0) ). 



(36) 



(37) 



R,o 



Notice that the last two current contributions, jA P ( r ') an d jAd( r ')> are written as a sum over 
augmentation sites and vanish outside the augmentation regions, where the all-electron and 
pseudo partial waves coincide. 

By construction, the current j^(r') computed within the GIPAW formalism is, as all 
physical observables should be, invariant upon translation of the system by a vector t, i.e. 
after translation the new current should be j^(r' — t). Interestingly, all three terms, jb2. e (r'), 
j^p(r') , and j^(r'), are individually invariant upon translation. The invariance of j^(r') is 
obvious from the definition of the AJ^(r') operator, Eq. (p7|) . The invariance of the other 
two contributions is less evident, and to prove it, we need to manipulate Eqs. (|35[) and (j36|). 
To this end, we notice that the second term in the r.h.s. of Eq. ( p5|) can be rewritten as a 
commutator, 



lpP s (r')B x r' = 2^i(^ )|i[B x r' ■ r, J p (r')] |*f>. 



(38) 



We can now use the generalized /-sum rule established in Appendix Eq. ( |A7|) , with the 
operators J p (r') and AJ^,(r') in the place of O and the operator r in the place of S, to 
rewrite the second terms in the r.h.s. of both Eqs. (|35|) and (1361), obtaining: 



( r > 

J bare V 



4^Re 



Bxr' 

2c 



•v|*<°>> 



, (39) 



(^ 0) |AJ R ,(r')^(. o )^ (1) |^ 0) ) - (^ 0) |AJ^(r')^(, o )5^ . v|^) , 

(40) 

where v = l/i[r, H^] is the velocity operator. Now the translational invariance of JbaL;( r/ ) 
and jAp( r ') is more explicit, since on translation both Bxr' ■ v/2c and B x R' ■ v/2c generate 
an extra term equal to B x t ■ v/2c, as does H^ l \ if, after the translation, we rewrite 
in terms of the variable of the translated coordinate system (r — t). 



jk 1 p(r')=4E^ 

R'.o 
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V. CURRENT RESPONSE IN EXTENDED SYSTEMS 



In Section [TIT] we developed a theory for a system containing a single augmentation 
region located at the origin, and then later for several augmentation regions. We must now 
check that our results are still useful in situations involving an infinite number of these 
augmentation regions, as is the case in the solid state. 

The expression for j^d(r') given by Eq. ( [371) can be straightforwardly applied to solid 
state calculations. But, the contributions to the all-electron current j^ TC {r') and jA P ( r ') 
given in Eq. (|39| ) and Eq. (|40| ) involve expectation values of the position operator. As these 
are not generally defined in an extended system, one might worry that Eqs. (|39|) and ([47)1 
are not valid. However, if they are rewritten in the following way: 



Jbarel 1 ) 



£Re 



(4f | J p (r')g(e ) (r - r') x p + £(R - r') x v£ ■ B|tf J») 



R 



(41) 



and 



R',o 



(^°)|AJ^(r')^( £o ) ((r - R') x p + £(R - R') x vg) ■ B|^) 



(42) 



then it becomes clear that they are indeed well defined. The Green Function operator G(s ), 
in an insulator, and both the paramagnetic augmentation operator AJ^,(r') and the non- 
local pseudopotential operator vg are short rang edE This ensures that contributions to the 
current response for large values of (r — r'), (R — r'), (r — R'), or (R — R') in Eqs. ([H]) and 
(f42T) vanish. 



VI. CURRENT RESPONSE IN INFINITELY PERIODIC SYSTEMS 

The expressions given above are valid for any extended system. However, the only such 
computationally tractable systems are those exhibiting translational symmetry, or infinitely 
periodic systems. We now develop the equations making this translational symmetry ex- 



plicit, by writing the electronic states as Bloch functions, |\E^°k) = e ik ' r |u^), where k is a 
reciprocal space vector within the first Brillouin zone and the corresponding eigenvalues are 
£ n ,\i- The cell-periodic function (rlu^jj is normalized within the unit cell. 

In order to take full advantage of this translational symmetry we first define the functions 
S b are(r',g) and S Ap (r', q) as: 



,ikT|„-;(0) 



'bare 



y, <?) 



E E Re 



i=x,y,z o 



R 



(43) 



q) 



E E Re 

i=x,y,z R',o 



(44) 
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where the Uj are unit vectors in the three Cartesian directions. We can then write 

jbarefr') = J™ ^ [Share (r', <?) ~ Sbaro(r', -q)] (45) 



j«(r') = Kmi[S Ap (r', g ) - S Ap (r', -q)] . (46) 

This can be seen to be correct by expanding the exponentials in Eqs. (|43 ) and fl44]) as 
e iq\ii-x = x _|_ . x + 0((gx) 2 ), taking the limits in Eqs. (f|5|) and ([5]) and comparing to 
Eqs. ([|l]) and (0). The limits taken using the expanded exponentials are valid since only 
finite values of x contribute to the total current (as established in Section |VJ) . 

The description of the electronic states as Bloch functions allows us to approximate the 
summations over the infinite number of occupied states in Section |V| as finite summations 
over k-dependent quantities. The k-dependent Green function is, 

li7 (0) Vi7 (O) l 

ft(e) = E ' _ ' ■ (47) 

A consequence of re-expressing the current contributions in terms of Sb a re( r/ , l) and SA p (r', q) 
is that we must evaluate several quantities at k and k + q simultaneously. For example, the 
usual form of the k-dependent non-local pseudopotential operator is generalized: 

C = EEIt)CCl' ( 48 ) 

r n,m 

This operator acts on Bloch functions at k to the left, and k' to the right. For k = k', 
V£y coincides with the k-dependent non-local pseudopotential operator, implemented in 
the plane- wave pseudopotential codes. The k-dependent projectors in terms of |pR, n ), the 
real space projectors, are given by 

fe) = E^ ik - (r ^ T) |PL + ^), (49) 

L 

where the L are lattice vectors and the r are the internal co-ordinates of the atoms. We 
arrive at analogous expressions for both the velocity operator, 

v k , k , = -*V + k' + -[r,K;U (50) 

and the paramagnetic current operator, 

jP k/(r0 = HV + k)|rO(r-| + |rO(r-|(-^V + kQ ^ ^ 

Combining the above we arrive at a compact expression for Sb are ( r/ , #) : 



>barc(r',g) = -4j- E E Re ^( U i°k|Jk,k+q I ( r ')^k+q I (£o,k)B X ■ V k+qi , k |w^) 

k i=x,y,z o,k 



(52) 
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where = guj and is the number of k-points included in the summation. Similarly, by 
also defining: 

AJ£, T , k , k ,(r') = £ [ft |B > [(0L + .,n|J p (r')|0L +T , m ) - (k+r,n\J p (r')\4> h+T , m )} (p£j, (53) 
the expression for the paramagnetic augmentation term is: 



2 



ciV k . — j — , 

K t=x,y,z L,r,o,k 



(54) 

These expressions for Sb a re( r ', <?) an d SA P (r', g) allow the evaluation of the all-electron current 
response through the Eqs. fl45|), fl46|) and 



VII. SUMMARY OF APPROACHES 

There are three different approaches that we could take in the calculation of the first 
order current response to a uniform external applied magnetic field. If the current response 
in an extended periodic system is required, then the approach described in Section |VI] must 
be taken. In this case, the expressions given in Eqs. P5j), (^BJ), ([52]), (|51|) and Q3"7| ) are 
evaluated, and it is referred to as the "crystal approach". The total current response in a 
finite system can be calculated using Eqs. (|35"D , (|3~6D, and (|3~7|). This approach is referred 
to as the "molecular approach" . Alternatively, using the results of the generalized /-sum 



rule, Eqs. (|39|), (fiO|), and (|37| ) can be used. This is the "molecular sum rule approach". 
Setting 7b = 1 in the GIPAW formalism, i.e. in the all-electron case, the "crystal approach" 
becomes equivalent to the MPL method,! the "molecular approach" becomes equivalent to 
the single gauge method (Eq. (3) of Ref. [H]), and the "molecular sum rule approach" 
becomes equivalent to the continuous set of gauge transformation methodH (CSGT) with 
the d(r) = r gauge function (Eq. (8) of Ref. [TT| ). 

The "crystal approach" can be used to calculate molecular properties through the use 
of large supercells. If the generalized /-sum rule holds, then the results obtained by each 
of the three approaches should be equivalent. This is demonstrated in Section [IX|. If the 
generalized /-sum rule does not hold well (for example, if the basis set used is far from 
completeness as is the case for the atomic-orbital basis sets used in most quantum chemical 
calculations^), then the "crystal approach" and the "molecular sum rule approach" will still 
give the same results. However, the results obtained using the "molecular approach" will be 
different. In particular, we expect that the "molecular approach" will require a much larger 
atomic-orbital basis set to converge the NMR chemical shifts than the other two methods, 
as it has been proved to be the case for all-electron Hamiltonians.il This is because the 
two terms in Eq. ( |35"1) (as well as the two terms in Eq. (|36|) ) of the "molecular approach" 
converge at different rates with respect to the completeness of the basis set .Hill 

VIII. CALCULATION OF NMR CHEMICAL SHIFTS 

It is important to show that the GIPAW method is a practical approach to the calculation 
of NMR chemical shifts. We have therefore implemented the method into a parallelized 
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plane-wave pseudopotential electronic structure code.l±3 Such codes self consistently calculate 
the ground state electronic structure. Specifically, the self consistent Hamiltonian and 
the corresponding wavefunctions \^^) that appear in the above expressions are obtained. 
In this section we outline the features of the implementation that are specific to the GIPAW 
method, and not to the pseudopotential method in general. The plane-wave pseudopotential 
method is most naturally suited to the "crystal approach" for the calculation of NMR 
chemical shifts. However, we also implemented both molecular methods in our plane-wave 



code, for completeness. This is described in Section VIII D 



A. Application of the Green function 

There are several points at which first order wavefunctions of the form, 

\^ ) )=G(e n )H {1) \^ ) ) (55) 
must be evaluated. The Green function Q(e) is given by, 

|t&(°)\/tjK°)| 

g(g) = E . (56) 

and a naive approach would require the explicit summation over all empty states. This is 
unnecessarily arduous. We can multiply Eq. ( |55D through by [e n — H^). If we then write 
Q — l^e°' ) )(^ , e ' ) l = 1 ~J2o I ^o^) (^o^ I > where the sums over o and e are over the occupied 
and empty states respectively, we obtain: 

(e„- J ffW)|*«) = eff«W> (57) 

This is a linear system involvingonly the occupied states, and can be solved using a conjugate 
gradient minimization scheme,^ as in Ref. [7[ This approach ensures that our method is 
comparable in computational cost to the calculation of the ground-state electronic structure. 



B. The velocity operator 

The velocity operator v = |[r, 5"^] appears in various guises throughout the relevant 
expressions above. The velocity operator may also be written as the first derivative of 
the k-dependent Hamiltonian with respect to k. The term related to the kinetic energy 
is straightforward to evaluate, and is simply the momentum operator. The term due to 
the non-local potential, which is defined numerically, is best obtained numerically. In our 
implementation we simply take the appropriate numerical derivative of the k-dependent 
non-local potential operators. The derivatives are evaluated by calculating the non-local 
potential at, say, k and k + q, where q is chosen to be small enough that the resulting 
numerical derivative is accurate, but not so small so as to introduce numerical noise. 
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C. The crystal approach 



The "crystal approach" requires that the limits of Eqs. (45) and ^ are evaluated. These 
are, in effect, similar to the numerical derivatives which must be take in reciprocal space 
in order to evaluate the velocity operators. And in practice we take the same value for the 
reciprocal space step size in both cases. The same considerations apply. The step should be 
chosen to be small enough that the resulting limit is accurately approximated, but not so 
small that numerical noise dominates. A typical value is 0.01 Bohr -1 . 



D. Finite systems in periodic boundary conditions 



Both the "molecular approach" and the "molecular sum rule approach" were imple- 
mented. The major difference between these approaches and the "crystal approach", from 
a computational perspective, is that the reciprocal space numerical derivative is replaced 
by a direct application of the position operator to the wavefunctions. Clearly, the position 
operator is not defined within periodic boundary conditions. But we can treat it approxi- 
mately by constructing a periodic saw-tooth like function (in practice we build the function 
in reciprocal space). Near the center of the simulation cell, or about wherever the saw-tooth 
is centered, this operator approximates the position operator. This approximation improves 
as the size of the simulation cell is increased, and for good results the magnitude of the 
induced current should be small on the surface where the saw-tooth function changes sign. 



E. From the current to the NMR chemical shifts 



The GIPAW approach separates the contributions to the current response into a bare 
term, j^eM; anc ^ two correction terms, the paramagnetic and diamagnetic corrections, 
j^p(r) and j^dM respectively To compute the NMR chemical shifts, using Eq. (Q), the 
induced magnetic field, B-£'(R), must be evaluated at each nuclear position R. In principle, 
one could combine the three current contributions and obtain (R) from the total current 
using Eq. ([]]). We use a different approach. We take advantage of the linearity of Eq. 
(HI), and we solve it for each of the three current contributions, obtaining a bare induced 
field B^ re (R), a paramagnetic correction field, B^ P (R), and a diamagnetic correction field, 
B^ d (R). 

To compute the correction fields, we suppose that just the correction currents, j^pM 

and jAdl 1 "); within the augmentation region f2 R contribute to B^(R) and B^ d (R) at the 
nuclear position R. Using this on-site approximation, combining Eqs. (|l]) and fl3?|), we 
obtain: 



B2(R) = 2 £ (^ 0) |p R , n )e^(p R , m |^(°)), (58) 

o,n,n 



where 



h tl . (R — r) x [B x (R — r)l . (R — r) x [B x (R — r)l ~ , , rn . 
< n = (Knl- ^iR-r 3 l0R,m) " ^ R ' n| 2c 2 lR-r 3 ( 59 ) 
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The coefficients n depend only on the atomic species, and need only be calculated once. 
Similarly, within the on-site approximation, by combining Eq. (P with the equations for 
the jAp( r ) correction current, we obtain expressions for the paramagnetic correction field, 
B^p(R), which depend linearly on the coefficients f^ n , 

f m,n = (0R,»l r L ^| 3 l0R,m) - (4>B.,n\ , 3 1 4>B.,m) , (60) 

where Lr = (r — R) x p is the angular momentum operator evaluated with respect to the 
atomic site R. Again, the coefficients f^ n depend only on the atomic species, and need only 
be evaluated once. 

To compute the bare induced field, B[ ) 1 a ' ) re (R), we Fourier transform Eq. (|l|) and j b a re ( r/ ) 
into reciprocal space. The induced magnetic field can then be simply evaluated as, 

B<1(G) = f ' GX f° (G) . (61) 

where G is a reciprocal lattice vector. We subsequently obtain B^^R) by a slow (since we 
only need the results at a few points in space) Fourier transform at the nuclear positions R. 
For G = 0, Eq. (^1|) can not be applied. Indeed the G = component of the induced 



magnetic field is not abnlk propertyB The G = component of the induced field is affected 
by the surface currents which appear on the surface of the sample. In particular, its value 
depends on the the shape of the sample, and is determined by macroscopic magnetostatics. 
Following the experimental convention, we assume a spherical sample in our calculations, 
for which: 

B i ( n ) (G = 0) = ^xB, (62) 

where x is the macroscopic magnetic susceptibility.0 To be consistent with the on-site ap- 
proximation for the correction currents, we should not take into account the contribution of 
j£ p (r) and j^J(r) to b£ } (G = 0), and so use: 

B«(G = 0) = ^Xba re B, (63) 

where Xbare is the contribution to the macroscopic susceptibility coming from the bare current 
j bare (r). Within the "crystal approach", we use the following ansatz for Xbar e : 

X^ J iq) - 2 " F{ °J + " F( - q \ (64, 
where Fij(q) = (2 — Sij)Qij(q), i and j are Cartesian indices, 

Q(v) = — E E Re [(^°k|u* x HV + k)£ k+qi (£ 0ik )ui x v k+qijk |^ k )] , (65) 

C k c i =X: y,z o,k 

and V c is the unit cell volume. In support of this ansatz, one can show that, when Tq = 1, 
i.e. in the all-electron case, the definition of Xbare, Eq. ( ]6l| ) becomes equal to the expression 
for the calculation of the all-electron macroscopic magnetic susceptibility, as derived in Ref. 
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F. Projectors 



In our implementation, we use norm-conserving Troullier-Martins pseudopotentialstd 
with single projectors for each angular momentum channel. As a result, the argument in 
Section [III B| holds and the bffl m terms are zero. However, in contrast to what Van de Walle 
and Blochl found for the calculation of hyperfine parameters,!^ we found that a minimum of 
two projectors per channel were required to ensure good transferability of the GIPAW cur- 
rent corrections. Otherwise, the projectors are constructed as described in Ref. [18], except 
that we choose a polynomial step function /(r) so that the pseudowavefunctions are cut off 
smoothly at some distance less than the pseudopotential core radius. 



IX. NUMERICAL TESTS OF THE GIPAW METHOD 

A. Comparison with IGAIM results 

Quantum chemical approaches have long been able to predict the NMR chemical shifts 
of small molecules, and one of the most widely used is the GAUSSIAN94fl quantum 
chemical code. Gregor et a/0 used this code to optimize the geometry and calculate the 
isotropic chemical shift of a selection of small molecules using both the GIAO and IGAIM 
methods.! We compare our GIPAW results (all chemical shifts reported here have been 
calculated within the Local Density Approximation^) to the IGAIM results for several of 
these molecules (using exactly the same relaxed geometries) in Table |I|. The total isotropic 
chemical shifts computed with GIPAW agree very well in all cases with the GAUSSIAN94 
results. 

The GIPAW results presented in Table | were evaluated using the "crystal approach" , 
but results obtained using the molecular approaches differ typically by less than 0.1 parts- 
per-million (ppm) in sufficiently large simulation cells, as demonstrated in Table 0. 

The GIPAW results are converged to the 0.1 ppm level using a plane- wave cut-off of 1QQ 
Rydbergs, a super-cell volume of 6000 Bohr 3 and a 2 x 2 x 2 Monkhorst-Pack k-point gridHI 
The states indicated in Table | were treated as core states in the pseudopotential calculations. 
The core contribution to the GIPAW chemical shifts is assumed to be constant (following the 
observations of Gregor et a/0), and evaluated in an all-electron atomic code. For hydrogen 
a pseudisation core radius of 1.2 Bohr was used and only the s-channel was augmented. As 
a result, since the paramagnetic correction term is proportional to the angular momentum 
of the augmentation channel (see Eq. ^), only the bare and diamagnetic correction terms 
contribute to the total isotropic chemical shifts. There is no core contribution for hydrogen. 
For the carbon shifts the s- and p-channels were augmented and a core radius of 1.6 Bohr 
used in the generation of the pseudopotential. For silicon and phosphorus the d-channel was 
also augmented and core radii of 2.0 Bohr used in both cases. Gregor et al attempted to 
converge the chemical shifts with respect to their localized basis set size, and the convergence 
appears to be to the 1 ppm level for the carbon and silicon shifts (see Fig. 2 of Ref. |T7D . 
However, the convergence appears to be less complete for the phosphorus shifts. It is just 
these chemical shifts for which the GIPAW and IGAIM results differ the most (although 
the errors as a fraction of the range of the chemical shifts are similar for all nuclei). While 
the diamagnetic correction term is found to be rigid with respect to chemical environment, 
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both the bare and paramagnetic correction terms are found to be strongly dependent on the 
system. The correction terms introduced by the GIPAW approach are therefore seen to be 
important even in the prediction of relative chemical shifts, and the rigid nature of the core 
contribution is reconfirmed. 

In Table |TTJ we examine the robustness of the GIPAW method with respect to pseudopo- 
tentials used. A variety of Troullier-Martins pseudopotentialsi^, with core radii ranging 
from 1.2 to 1.8 Bohr, were used to calculate the NMR chemical shift for carbon in methane. 
While the bare contribution to the chemical shift is observed to change by over 10 ppm, the 
total shifts, including the GIPAW correction terms, are constant to within 1 ppm. There is 
virtually no difference in the total shifts between potentials with core radii of 1.2 and 1.4 
Bohr. 



B. Comparison with all-electron plane- wave results for diamond 

As the GIPAW method presented here is, to the authors' knowledge, the only approach 
available for the calculation of all-electron NMR chemical shifts in solids, a truly independent 
validation is not possible. However, by constructing a suitable pseudopotential and taking 
a high enough plane-wave cut-off energy we are able to compare with essentially all-electron 
results — in which all the electrons in the chosen system are considered to be valence 
electrons. In this way we can check the corrections to the conventional pseudopotential 
results. Obviously, such calculations are computationally intensive due to the extremely 
large number of plane-waves required to reach convergence. We therefore choose diamond 
as our example periodic system. Carbon is sufficiently light that an all-electron plane-wave 
calculation is possible, and the diamond structure has a very small primitive unit cell and 
a high degree of symmetry. The Is, 2s, and 2p electrons are all considered to be valence 
electrons and we construct a purely local Troullier-MartinB pseudopotential with a core 
radius of 0.4 Bohr radii. 

Table [TV] compares the results of a GIPAW pseudopotential calculation (the Is electrons 
are treated as core electrons, and a core radius of 1.6 Bohr radii used) and the all-electron 
plane-wave calculation obtained with the purely local Troullier-Martin pseudopotential. The 
contributions can be separated into core and valence terms in a gauge invariant way, as shown 



in Ref. [IT]. Thus, in the case of the all-electron result we performed two calculations of the 
chemical shift after achieving self-consistency, once taking into account all the electrons, 
and a second time excluding the valence electrons from the calculation of the chemical shift. 
The valence term presented is the difference between these two results. We present the 
all-electron results at two plane-wave cut-offs — 800 and 1400 Rydbergs and a 10 x 10 x 10 
Monkhorst-Pack k-point grid. All the contributions to the chemical shifts are converged to 
within a part-per-million. The valence contributions of the GIPAW and all-electron results 
differ by only 1.39 ppm which may be attributed to the slight uncorrected pseudisation 
error that remains in the all-electron result. We have confidence that if the core radius were 
reduced to less than 0.4 Bohr radii the difference between the results of the two approaches 
would decrease. The GIPAW pseudopotential result is expected to be closer to the true 
all-electron NMR chemical shift. 
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X. CONCLUSIONS 



We have presented an ab initio theory for the evaluation of NMR chemical shifts in both 
finite and infinitely periodic systems. We have correctly treated the complications introduced 
due to the use of pseudopotentials, and so, in contrast to the original implementation of the 
MPL approach,! we are not restricted to the calculation of the chemical shifts for light 
elements. We introduced an extension to the Projector Augmented- Wave method which 
is valid for systems in non-zero uniform magnetic fields, the Gauge Including Projector 
Augmented- Wave method. 

Our implementation of GIPAW into a parallelized plane-wave pseudopotential code al- 
lows the calculation of NMR chemical shifts in large, low symmetry extended systems. We 
expect that the methodology will prove useful in the calculation of other magnetic prop- 
erties. Our work also suggests that the implementation of GIPAW into quantum chemical 
approaches would lead to a considerable improvement in their efficiency for the calculation 
of NMR chemical shifts for heavy elements. 
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APPENDIX A: THE GENERALIZED F-SUM RULE 



The generalized /-sum rule holds for any pair of hermitian operators O and £, where O 
and £ are respectively odd and even on time reversal, i.e.: 



and 



\ow) = -W\o\<f) 



{<j>\£W) = ms\ 



(Al) 



(A2) 



for any \<p) and \<f>') such that (r\<p) and (r|0') are real. It is straightforward to verify that 
p, L, v, Vj£, J p (r'), and AJ^(r') are odd, and that r and operators that are a function of r 
are even. To derive the sum rule, we consider the quantity 



8 = -4$>e {yW\og(e )-[£,HM]\V) 



(A3) 



The sums over o and d (below) run over the occupied orbitals, and those over e' over 
the empty ones. Using the fact that H^\¥^) = e k \¥^), Eq. (|33]) and £ e , |^ e ')(^e'l = 
1 — So' l^o'X^o'l the expression for s may be rewritten as, 



-4^Re +4^Re[i(f(°)|O|fi? ) )(fJ ) |£:|f( )) 



(A4) 
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Since the eigenstates l^j^ ) can be chosen in such a way that (rl^^) is a real quantity, 



(0)\ 



(¥^\0\^) = -(#H0|*n and i^k'l^k' 1 ) = WVl^D- Using these relations it 
follows that: 



,(o)\ 



,(o)\ 



,(0)|ChT,(0)\ 



r(0)|PliT,(°)\ 



0,0' 0,0' 

= -Y,(mo\&°WS ) \e\W)> 



(A5) 



0,0' 



where for the last equality we just interchanged the dummy indexes o and 0' . From Eq. 
we conclude that the double summation of Eq. (|A4|) is equal to zero and: 



(A6) 



a = -4]TRe f-(^ 0) |^|^ 0) )l =2j2(^ ) \-[£,0]\¥^). 



^From this expression we finally obtain the generalized /-sum rule: 

2E(^ 0) I-[^]I^ 0) > = ~ A J2 Re bi 0) l^(^o)-[^,^ (0) ]|^ 0) ) 

o 1 n L 1 



(A7) 
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TABLES 



TABLE I. Isotropic absolute chemical shifts calculated using the IGAIM method by Gregor et 
and the corresponding GIPAW-LDA results. The GIPAW calculations were performed using 
a plane- wave cut-off of 100 Ry and in a 6000 Bohr 3 simulation cell. With "bare", "Ad", and 
"Ap" , we indicate the valence GIPAW contributions to the chemical shifts, given by the bare field 
B^2. e (R) and the two correction fields B^(R) and B^p(R), respectively. The core contribution 
to the GIPAW chemical shifts is assumed to be constant and evaluated in an all-electron atomic 
code. All quantities are given as ppm. 



Molecule 


Core 


bare 


Ad 


Ad 


Total 


Total 


T~T atnin 

11 CJj U VJ11J. 














CH 4 


0.00 


30.47 


0.40 


0.00 


30.87 


30.99 




0.00 


25.71 


0.41 


0.00 


26.13 


26.50 


K jd H ft 


0.00 


22.33 


0.41 


0.00 


22.74 


23.25 


TMS 


0.00 


30.41 


0.40 


0.00 


30.80 


31.02 


SiH 3 F 


0.00 


24.92 


0.38 


0.00 


25.30 


25.13 


SI2H4 


0.00 


24.53 


0.36 


0.00 


24.90 


24.78 


SiH 4 


0.00 


26.96 


0.37 


0.00 


27.33 


27.28 


O atom 


Is 












CO 


198.88 


-126.25 


4.59 


-100.15 


-22.93 


-21.16 


CH 4 


198.88 


16.86 


3.97 


-28.76 


190.96 


191.22 


CHcF 


198.88 


-49.64 


3.93 


-54.70 


98.47 


99.66 




1 98 88 


-1 3 98 


3 91 


-39 05 


1 49 77 


1 50 44 




198.88 


-89.51 


4.07 


-77.32 


36.12 


39.52 


CF 4 


198.88 


-92.12 


3.51 


-76.05 


34.22 


35.29 


TMS 


198.88 


9.12 


3.97 


-32.65 


179.33 


182.08 


Si atom 


Is2s2p 












SiF 4 


832.39 


-19.43 


5.28 


-408.26 


409.97 


409.69 


SiH 3 F 


832.39 


-19.50 


5.70 


-510.30 


308.29 


305.45 


Si2H 4 


832.39 


-9.04 


5.80 


-622.45 


206.70 


202.99 


SiH 4 


832.39 


-0.21 


5.98 


-410.20 


427.97 


424.37 


TMS 


832.39 


-17.39 


5.70 


-518.00 


302.70 


304.39 


P atom 


Is2s2p 












PF 3 


902.47 


-32.94 


6.08 


-697.61 


178.00 


172.52 


P 2 


902.47 


-33.84 


7.58 


-1236.95 


-360.75 


-375.45 


P 4 


902.47 


49.84 


7.42 


-126.79 


832.94 


826.62 
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TABLE II. Comparison of the three different GIPAW approaches described in Section VD. 
The GIPAW-LDA calculations were performed using a plane- wave cut-off of 100 Ry and in a 6000 
Bohr 3 simulation cell. The total isotropic chemical shifts are given as ppm. 



Molecule 




Molecular 




Molecular sum rule 




Crystal 


H atom 
















CH 4 




30.75 






30.76 




30.87 


CH 3 F 




26.02 






26.01 




26.13 






22.69 






22.69 




22.74 


TMS 




30.76 






30.76 




30.80 


S1H3F 




25.40 






25.40 




25.30 


Si 2 H 4 




24.92 






24.93 




24.90 


SiH 4 




27.57 






27.58 




27.33 


C atom 
















CO 




-22.92 






-22.90 




-22.93 


CH 4 




191.08 






191.09 




190.96 


CH 3 F 




98.53 






98.52 




98.47 


CH3NH2 




149.61 






149.62 




149.77 






36.13 






36.14 




36.12 


CF 4 




34.62 






34.30 




34.22 


TMS 




179.17 






179.19 




179.33 


Si atom 
















SiF 4 




410.12 






409.85 




409.97 


SiH 3 F 




308.27 






308.23 




308.29 


Si2H 4 




206.50 






206.49 




206.70 


SiH 4 




427.95 






427.95 




427.97 


TMS 




302.61 






302.61 




302.70 


P atom 
















PF 3 




177.90 






177.70 




178.00 


P 2 




-360.97 






■360.97 




-360.75 


P 4 




832.87 






832.87 




832.94 


TABLE III. 


The NMR chemical shift for carbon 


in methane using Troullier-Martins potentials 


with a range of 


core radii. These LDA calculations 


were performed using 


a plane-wave cut-off of 


180 Ry (convert 


;ed to 0.01 ppm for the hardest potential) and in a simulation cell of 1000 Bohr 3 . 


With "bare", "Ad", 


and "Ap" , we indicate the valence GIPAW contributions to the chemical shifts, 


given by the bare field B| )are (R) and the two correction fields B^ 


(R) and Bg(R), 


respectively. 


Core radius (Bohr) 






CGIPAW 












Core 


Bare 


Ad 




Ap 


Total 


1.2 




198.88 


7.30 


3.96 




-19.14 


191.00 


1.4 




198.88 


12.22 


3.99 




-24.08 


191.01 


1.6 




198.88 


17.03 


3.98 




-28.64 


191.25 


1.8 




198.88 


21.65 


3.92 




-32.86 


191.59 
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TABLE IV. The valence contribution to the isotropic chemical shift of crystalline diamond 



(ppm). 



Method 


Valence contribution to a 


GIPAW 


-65.85 


All-electron at 800 Ry 


-64.89 


All-electron at 1400 Ry 


-64.46 
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